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PACS. 62.20.Qp - Tribology and hardness. 

PACS. 68.35.Ct - Interface structure and roughness. 

PACS. 91.60.Ba - Elasticity, fracture and flow. 
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Abstract. - We model numerically the partial normal contact of two elastic rough surfaces 
with highly correlated asperities. Facing surfaces are unmated and described as self-affine with 
a Hurst exponent H . The numerical algorithm is based on Fourier acceleration and allows 
efficient simulation of very large systems. The force, F, versus contact area, A, characteristics 
follows the law F ~ j\( 1 + H >/ 2 m accordance with the suggestion of Roux et al. (Europhys. Lett. 
23, 277 (1993)). However finite size corrections are very large even for 512 x 512 systems where 
the effective exponent is still 20% larger than its asymptotic value. 
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The mechanical properties of bodies in contact have been studied for a long time |lj since 
they have important applications ranging from the flow properties of powders to earthquake 
dynamics J2|]. For example, frictional properties are related to the normal stresses that develop 
during contact, and they have recently been shown to be very sensitive to heterogeneities in 
the normal stresses Q . On the fault scale, the dynamical stress field, which is responsible for 
earthquakes, is strongly influenced by heterogeneities due to asperity squeeze [|l|. 

In this letter we investigate numerically the elastic response of self-affine rough surfaces 
that are squeezed together. The rough surface, which we take to be oriented in the (x, y) 
plane, is given by z = z(x,y). If p(z;x, y) is the probability density to find the surface at a 
height z at (x,y), self affinity is the scaling property 

X H p(X H z; Xx, Xy) = p(z; x, y) , (1) 

where H is the Hurst or roughness exponent (q,H|. There is strong experimental evidence that 

surfaces resulting from brittle fracture are self affine with a Hurst exponent equal to 0.8 J7|-p^|. 

When two elastic media are forced into contact along two non-matching self-affine rough 

surfaces, two mechanisms conspire to produce a power law dependence of the applied normal 
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force F on the penetration depth, D: (1) the contact area increases as D is increased, and 
(2) the fluctuations in the normal stress field where there is contact reflect the self affinity of 
surfaces. In the much simpler case of a spherical asperity in contact, Hertz showed about 120 
years ago |L3[ that F ~ Z? 3 / 2 , and A ~ D. The Hertz law has been verified experimentally 
in the case of a diamond stylus sliding on a diamond surface E3 . The power law dependence 
of the force on penetration has also been established experimentally for fracture surfaces in 
contact ]l5| ]. 

In 1993, Roux et al. p6[ proposed the following dependence of F and A on D for self-affine 
elastic surfaces in contact, 

F ~ D 1+1 ' H , (2) 

and 

A ~ D 2 / H , (3) 

where H is the smallest Hurst exponent of the two surfaces. The argument goes as follows. 
First, we note that the rough surface with the smallest Hurst exponent is the roughest and 
will dominate the scaling properties of the common interface of the two media in contact up 
to a cross-over length scale |17|| . Second, we note that nowhere in the arguments that follow 
will we need both surfaces to deform. Thus, we will make the assumption that the surface 
with the largest Hurst exponent (or one of them if they have the same Hurst exponent) simply 
is flat and elastic, whereas the other one is rough and infinitely rigid. Let us now rescale the 
spatial coordinates of the surface, x — ► Xx, y — > Ay and z — > A z. As a consequence of Eq. 
(pi), the surface is statistically invariant under this operation. The penetration D needs to be 
rescaled D — » X H D as it "points" in the z direction, and the contact area A — > X 2 A as it lies 
in the (x,y) plane. The local deformation, u, of the surface must be rescaled as u — > X H u, 
since u also "points" in the z direction. The local deformation, u, at (x,y) is related to the 
normal stress field a by the expression |l3| 

u(x,y)= dxdyG(x-x',y-y')a(x',y') , (4) 

where G ~ 1/r, and r = \(x — x',y — y')\. Thus, we find that a — » X H ~ 1 a under rescaling |18| . 
The total force F — A<r, and consequently, F — > X 1+H F under rescaling. Eliminating A 
between any pair of variables in the above expressions, gives the power law dependence of the 
different variables. In particular, we find Eqs. (||) and (||) and 

F ~ A( 1+ff )/ 2 . (5) 

Eq. (g) was tested experimentally and numerically in Ref . |1^] . The obtained results were 
consistent with the suggestions of Roux et al. (TsJ , but hardly convincing. The first problem 
was that the definition of the zero level of the penetration is not the point of first contact 
between the two surfaces, but rather has to be treated as a free parameter in the data fits. 
Another numerical problem encountered was the necessity to invert dense matrices, with the 
consequence that only small systems could be studied (up to 32 x 32). 

In the present letter, we have developed a new and much more efficient algorithm that 
allows us to simulate easily systems up to 512 x 512. We are, therefore, able to present 
convincing tests of the scaling relation Eq. (|^). Whereas the numerics of Ref. |19| were based 



on six samples of size 32 x 32, we base our data on 300 samples of size 512 x 512. One surprising 
result in the present study is how strong the finite-size corrections to the results are. Size is 
estimated in units of the lower cut-off scale for the self-affine invariance of the rough surfaces. 
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Fig. 1 - A stiff medium with a rough surface (top) is pushed into an eiastic, fiat medium (bottom). 
The dashed, horizontai fine represents the undeformed medium at a given penetration depth, whiie 
the fuli fine represents the actuai surface of the medium. The penetration depth D is defined with 
respect to the dashed line. 



We now describe the model and its numerical solution. We represent both surfaces on 
two-dimensional L x L lattices with deformations taking place only in the third transverse 
direction. The discretization size is adjusted to the unit size that is the lower cut-off length 
of the self-affine scaling. As the hard rough surface is pushed into the elastic one, the forces 
and deformations are described by the discrete version of Eq. m) , 



= E G ^< ( 6 ) 



with the Green function given by 



G,„ = ^ft°K (7) 

Zire \i — j\ 

Ui is the deformation of the elastic body at site i and fi the force acting at that point. In the 
Green function, s is the Poisson ratio, e the elastic constant, and \i — j\ the distance between 
sites i and j. The indices i and j run over all L 2 sites. 

To define the problem completely, the boundary conditions need to be specified. Clearly, 
in the regions where the rough surface is in contact with the elastic one, the deformation, Ui, 
is specified (it conforms to the shape of the rough surface) and one solves for the force. On the 
other hand, in the regions where contact has not been established, the elastic surface deforms 
in response to the influences from the contact regions. In this case equilibrium is established 
when the net forces acting on the surface vanish. Therefore, in the no-contact region, the 
force is specified, it vanishes, and one solves for the deformation. 

We can build these boundary conditions into the equations to facilitate solving them. To 
do this, we define the diagonal L 2 x L 2 matrix, K, with elements equal to 1 on contact sites 
and on free (no-contact) sites. Clearly the vector Km is zero everywhere there is no contact. 
At contact points, Ku is equal to the imposed deformation given by the shape of the rough 
surface. 

Eq. @ can be rewritten as, 

G(I-K)/+GK/ = (I-K)u + Ku, (8) 

where we use matrix-vector notation, and I is the identity. This form is convenient because 
as mentioned above, Km is a known quantity (boundary condition). In addition, the vector 
(I — K)/ is always zero because the force / is nonzero only at contact points. Putting the 
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unknowns on the left hand side and the boundary conditions on the right hand side of Eq. (0) 
we obtain, 

GK/ (I - K)u = Kw . (9) 

Now define the vector x representing all the unknown quantities. Clearly, 

x = K/+ (I - K){? . (10) 

With this definition, and noting that K(I - K) = (I - K)K = 0, K 2 = K, and (I - K) 2 = 
(I — K) we can write Eq. (g) as 

GKx-(I-K)x = Ku, (11) 

and finally 

(I-(I-G)K)# = K«. (12) 

Eq. (|l2j) is of the familiar form, Af = b which can be solved using, for example, the Conju- 
gate Gradient (CG) method |2^,|l[- The difficulty is that the Green function G is represented 
by a full L 2 x L 2 matrix. This will lead to the number of operations per iteration scaling 
as L which is prohibitive. We overcome this difficulty by doing the matrix multiplications 
involving G in Fourier space by using FFTs. This leads to a stable and extremely efficient 
algorithm, scaling like L 2 ln(L), with which we were able to study systems with size up to 
L= 512. 

One more technical detail remains. The matrix K, which indicates the contact points 
needs to be determined. The problem is that as we push the rough surface into the elastic 
one, the latter deforms. Therefore, the contact area is not equal to the area obtained by 
simply taking a cut through the rough surface. We obtain the correct contact area as follows. 
Our initial assumption is that the contact area is equal to the area obtained from a simple 
cut of the rough terrain, see Fig. fy. This determines the initial K which is then used in Eq. 
(O) . The solution thus obtained gives the forces where there is contact and the deformations 
where there is none. Some of the forces thus obtained are negative since the elastic surface is 
trying to pull away from the rough surface. We therefore modify K by zeroing the elements 
corresponding to sites where the force is negative, and we solve again. We repeat this process 
until there are no sites with negative forces. This algorithm always converges giving the 
correct contact area and forces. 

To verify the algorithm and program we first tested it with a Hertz contact. Figure |2| 
shows the force-contact area characteristics of a hard sphere with radius R = 6000 which is 
pushed into an elastic, flat medium with elastic constant e = 100. The calculation was done 
on a 512 x 512 lattice. The exact Hertz solution gives A ~ F 2 / 3 , while a least-squares fit on 
the data gives A ~ p - 63 The exponent given by our model, 0.63, is in excellent agreement 
with the exact value, the difference being due to finite size effects. We have verified that for 
smaller systems and smaller asperities, the exponent moves farther away from the 2/3 value. 

Having established that the algorithm works correctly, we then simulated various aspects 
of squeezing a stiff self-affine surface into an initially flat elastic one. We measured the total 
force, F, as a function of the total contact area, A. In addition, for all situations simulated, 
we measured the force at each contact point thus obtaining a very detailed picture of the 
squeezing process. 

We show in Fig. || the force, F, versus contact area, A : for different system sizes and Hurst 
exponent H = 0.6, 0.8. The data are averaged over 300 samples for each size and H . The 
elastic constant for the elastic, flat medium is e = 100. Even though the force versus contact 
area curves produce very high-quality fits to power laws, F ~ A a , the resulting exponents, a, 
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Fig. 2 - Contact area versus force for a hard sphere with radius R — 6000 squeezed into an elastic flat 
medium with elastic constant e = 100 (i.e. a Hertz contact). The calculation is done on a 512 x 512 
lattice. The least-squares fit, assuming a power law gives A ~ p°- 63 



depend strongly on the system size, L, even for systems as large as 512 x 512. For H = 0.8, 
we find from Fig. |, a(L = 32) = 1.33, a (64) = 1.27, a(128) = 1.21, a(256) = 1.16, and 
a(512) = 1.12. Correspondingly, from Fig. y, we find for H = 0.6, a(L — 32) = 1.31, 
a(64) = 1.24, a(128) = 1.19, and a(256) = 1.13. 

In figure we show these exponents versus 1/ L b for H = 0.8 (circles) and H — 0.6 
(triangles). It is clear from this figure that the finite size effects are large. It is not easy 
to extrapolate to L — > oo using a(L) = a^ + aL~ b because the power law has a small 
exponent, b, and requires much larger sizes to give reasonable results to this three parameter 
fit. Instead, we will verify the result of reference [|16| by fitting to a one parameter function, 
a(L) = (1 + H)/2 + L~ b . The results are shown as solid lines in Fig. [|. The high quality of 
the fits strongly supports the result ctoo = (1 + H)/2. We found b = 0.24 for H = 0.8 and 
b = 0.2 for H = 0.6. It is not clear if the difference between the two exponents is meaningful. 
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Fig. 3 - left: Force, F, versus contact area, A, for 300 independent samples of hard self-affine surfaces 
with H — 0.8 pushed into an initially flat medium with elastic constant e = 100. System sizes from 
bottom to top: 32, 64, 128, 256 and 512. right: Same but with H — 0.6 and sizes up to 256. 
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Fig. 4 - Effective exponent a versus 1/L b , where L is the system size. Circles are for H — 0.8 and 
give 6 = 0.24, triangles are H — 0.6 and give 6 = 0.2. Curves are given by a — (1 + H)/2 + L~ b . 

It is worth noting that although the case of the Hertz contact did display some finite size 
effects, they were very small especially when compared with the large effects observed with self 
affine surfaces. This can be attributed to the long distance correlations in self affine surfaces. 

Thanks to a fast algorithm based on Fourier acceleration we were able to describe precisely 
the evolution of the contact between two elastic rough surfaces from the first contact of 
one asperity to the maximum contact (see below). Asperity roughness is characterized by 
a self-affine (Hurst) exponent in a range compatible with fractured surfaces. Our algorithm 
easily allows extension to other roughness exponents. We demonstrated that the theoretical 
prediction from Roux et al |lq| is valid but very strong finite size effects exist. For instance, 
the effective relationship between normal force and contact area is 



F ex ^ 1+ff )/ 2+z 



(13) 



were b rj 0.2. This finite size effect results in a large sensitivity to the cut-off scale of the 
self- affine scaling. 

Finally, we remark that when we push a hard rough body into the elastic initially flat one, 
the maximum contact area, A/L 2 , obtained is never unity! In Fig. |we show the total force, 
F, versus the contact area, A/L 2 , up to the maximum attainable contact area. It is clear from 
the figure that it is not unity: For the L = 512 system it is of the order of 20%. If we allow 
plastic yield, full contact would be attainable. This is straightforward to do in our model. 
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